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' Abstract. 

^ . We investigate the behavior of the moments of the partially transposed reduced 

^ I density matrix in critical quantum spin chains. Given subsystem A as union 

of two blocks, this is the (matrix) transposed of pA with respect to the degrees of 
, freedom of one of the two. This is also the main ingredient for constructing the 

C/2 ' logarithmic negativity. Wc provide a new numerical scheme for calculating efficiently 

all the moments of using classical Monte Carlo simulations. In particular we study 
several combinations of the moments which are scale invariant at a critical point. Their 
behavior is fully characterized in both the critical Ising and the anisotropic Heisenberg 
' O . XXZ chains. For two adjacent blocks wc find, in both models, full agreement with 

^ I recent CFT calculations. For disjoint ones, in the Ising chain finite size corrections 

^ ' are non negligible. We demonstrate that their exponent is the same governing the 

unusual scaling corrections of the mutual information between the two blocks. Monte 
Carlo data fully match the theoretical CFT prediction only in the asymptotic limit of 
^ I infinite intervals. Oppositely, in the Heisenberg chain scaling corrections are smaller 

and, already at finite (moderately large) block sizes, Monte Carlo data are in excellent 
agreement with the asymptotic CFT result. 
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1. Introduction 

In recent years there has been a growing interest in characterizing the behavior of 
entanglement related quantities (see [1] for general reviews) in many body quantum 
systems. Moreover, the deep connection between entanglement and conformal field 
theory (CFT) [1, 2, 3, 4] has boosted a huge amount of work at the frontiers between 
quantum information, condensed matter, and quantum field theory. 
Entanglement related quantities are usually constructed considering a bipartition of a 
system S into two subsystems as S = AU B. Given the (pure) state of the total 
system and the density matrix p = the reduced density matrix for subsystem A 

is obtained by tracing over the degrees of freedom of B a.s pa = Tr^ p. The entanglement 
between A and B can be quantified using the von Neumann entropy S'yi = — Tr pA log pa- 
Alternatively, from the n— th moment Tr(p^)(n G N) of the reduced density matrix one 
can construct the so called Renyi entropies S*^^ = l/{n — 1) logTr(p^)", which are 
also standard entanglement measures. The von Neumann entropy is recovered from the 
Renyi entropies as the analytic continuation Sa = lini„^.i S^^\ 

Let us now consider a ID critical quantum system (spin chain) described in the 
scaling limit by a conformal field theory (CFT). After restricting to periodic boundary 
conditions, and taking subsystem A a single interval (Fig. 1 (a)), the scaling behavior 
of the Renyi entropies is given as [2, 3, 4, 5, 6] 



with i the length of the interval. Here c is the celebrated central charge [7], a an 
ultraviolet cutoff (lattice spacing), and a non universal constant. From (1) the von 
Neumann entropy Sa is obtained as 



with Ci also non universal. 

Both (1)(2) are nowadays accepted as the standard tools to extract the central charge 
in ID critical systems, while other universal featmes can be obtained by analyzing their 
finite size corrections [8, 9, 10, 11, 12, 13]. 

From the field theory point of view, much more information about the underlying CFT 
is contained in the mutual information between two intervals [14, 15, 16, 17, 18, 19, 20, 
21, 22, 23, 24, 25, 26, 27, 28]. In fact, given A as sum of two non complementary blocks 
as A = AiU A2 (Fig. 1 (c)), their mutual information /^"^.^^ = 5*^"^ + S"^"^ — 5'^"u^2 
gives access to the full operator content of a CFT [29]. From the quantum information 
perspective, however, since subsystem A is not in general in a pure state, Iai-A2 "'^^ ^ 
measure of the mutual entanglement between Ai,A2, although it contains information 
about the correlations between them. A standard measure of the entanglement between 
two blocks is instead the so-called logarithmic negativity [30]. Denoting a basis for the 




(1) 



S'yi = - log - + ci 



(2) 



a 
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Figure 1. Different ways of partitioning a ID spin chain (of length L) considered 
in this work. Subsystem A is denoted with the thicker hne. (a) A is the single 
interval A = [ui, vi] (of length £ ^ \vi — ui\). (b) A is made of two adjacent intervals 
A = = [wi, t'i]U[u2, 1)2] (of length respectively i!i = wi| and f 2 = 1^2 — M2|). 

(c) Two disjoint intervals Ai,A2 at distance d = \u2 — Vi\. In this work we always 
considered intervals of equal length, i.e. ii — £2 — £■ 



Hilbert space of Ai{A2) as le^-^^) {\ef^)), one first defines the partially transposed reduced 
density matrix (with respect to A2) as 

(ef^ ® ef |p-^|e« ® ef )) = (ef ^ ^ ePlp^e^^ ® ef ) (3) 

Then the logarithmic negativity is readily given by 

^^^logllp^^lli (4) 

with IIpJIi denoting the trace norm, i.e. the sum of the absolute values of the 
eigenvalues of . Besides its importance in quantum information, it has been pointed 
out recently that the logarithmic negativity is a universal quantity at a second order 
phase transition [31], which makes S a usable indicator of critical behavior in quantum 
many body systems. 

Arguably in the context of field theories a major role is played by the moments of 
pj, i.e. Tr(pJ)". In particular, the recent CFT approach developed in Ref. [32, 33] 
provided full analytical understanding of their scaling behavior, for both adjacent and 
disjoint intervals (Fig. 1 (b) and (c)). In the former case this allowed, via the analytic 
continuation n — )■ 1, to obtain an exact expression for S. For two disjoint intervals, 
although it was not possible to perform the analytic continuation, it has been proven 
rigorously that at a second order phase transition £^ is a universal function of the 
(dimensionless) harmonic ratio 

~ {U2 - Ui){v2 - Vi) 

with Ui^Vi the endpoints of the two blocks (cf. Fig 1). Yet, we should mention 
that the analytical treatment of is a hard task (results are available only for free 
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bosons [34] [32, 33]), and, as a matter of fact, a precise verification of the above 
mentioned CFT findings in microscopic models was up to now lacking. 
In this work we demonstrate that all the moments of can be efficiently computed in 
classical Monte Carlo simulations exploiting the mapping between a quantum system 
in d dimensions and a classical one in d + 1. The Monte Carlo technique we propose, 
which is itself one of our main results, generalizes the one used in Ref. [15, 21, 22, 35] 
to compute Trp^ (and the Renyi mutual information thereof), allowing us to provide a 
robust verification of the CFT results in Ref. [32, 33]. To this purpose, following [32, 33], 
it is convenient to define for adjacent blocks the ratio r„ | as 



T2=L/4^„ 



IUA2) 



(6) 



Tr(PATuA2 

with z = i/L. Here the notation Trp^^~^^^^^^ means that the partial transposition is 
done with respect to the degrees of freedom of block A2 of length i{L/4). For two 
disjoint blocks we use instead the ratio Rn{y) 

- (7) 

PA1UA2 

Remarkably both r„ and Rn are scale invariant quantities at a second order transition 
(cf. [32, 33] or section 2), which makes them good indicators for quantum and classical 
critical behaviors. In particular, a part from scaling corrections, which are expected 
to be less severe for S [32], they are as effective as the logarithmic negativity. Also, 
their being not related to any local order parameter makes them suitable especially 
for detecting topological transitions. On the CFT side, we anticipate here that (cf. 
section 2), while the scaling function r„(2) is fully characterized in terms of the central 
charge, Rn{y) is a universal function of y and depends on the full operator content of 
the given theory. In this sense Rn{y) provides yet another tool (besides the mutual 
information) to unveil the deep structure of CFTs. 

The models. In this work we focus on the critical Ising quantum spin chain and the 
anisotropic Heisenberg XXZ model at A = —l/\/2 (A is the anisotropy). In both cases 
we consider periodic boundary conditions. The Ising chain in a transverse field h is 
defined in terms of the Hamiltonian 

'^' = -YJ^t<y-+i + h<TV, (8) 

i 

with erf '^'^ the Pauli matrices and i = 1,2, The model exhibits a ferromagnetic 
(paramagnetic) phase for h < 1 [h > 1) with a second order phase transition at the 
critical value = 1. Its critical behavior is described in the continuum by the free 
Majorana fermion theory, which is the simplest and most studied CFT (with c = 1/2). 
The same theory describes the critical behavior of the 2D classical Ising model. 

X Note that, at difference with Ref. [32], r„ is defined here without the logarithm. 
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The anisotropic Heisenberg XXZ spin chain is instead defined by the interaction 

L L 

H''^' = af^, + af+J + A ^ af^, (9) 

4 = 1 1 = 1 

Its phase diagram shows a critical hquid phase for — 1 < A < 1 and a gapped one 
at |A| > 1 (the point A = — 1 is critical but not conformal invariant). The liquid 
phase is described in the continuum limit by the free compactified boson theory (or 
Luttinger liquid), which is a conformal field theory with c = 1 [36, 37]. The region 
at — 1 < A < — 1/a/2 is also mapped into the low-temperature critical phase of the 
2D classical XY model, which describes a system of interacting classical spins (rotors) 
Si = (cos 6*4, sin 6*4) and is defined by the Hamiltonian 

(ii) (ij) {ij) 

Here ipi = e*^' G f/(l) are phases living on a two dimensional square lattice, /3 = 1/T, 
and (ij) denotes nearest-neighbor sites. The XY model exhibits a low-temperature 
gapless critical phase characterized by quasi-long-range order (QLRO). This is divided 
from the standard paramagnetic phase at high temperature by a Berezinskii-Kosterlitz- 
Thouless (BKT) topological transition at Pbkt = 1.1199(1) [38, 39, 40, 41, 42, 43]. The 
critical properties at the BKT point are the same as in the XXZ chain at A = -1/^2, 
a part from logarithmic corrections that are present only in the classical model. 

Summary of the results. The main results of this work can be summarized as follows. 
For two adjacent blocks, in both the critical Ising and XXZ chains and already for finite 
(large) i, the ratio rn{z) is numerically indistinguishable from its asymptotic value (i.e. 
at L,£ — J- 00), meaning that scaling corrections are small. Moreover, Monte Carlo data 
are in full agreement (for any value of z) with the CFT result in Ref. [32]. 
For disjoint intervals we focus on Rsiy). For the Ising chain unusual (in the 
sense of Ref. [10]) scaling corrections are non negligible, as observed for the mutual 
information [17, 21, 22, 23, 27] (see also [32, 33]). We numerically demonstrate that 
they decay as with co^ = 1/3, in agreement with the general behavior (as £~^/") 
found in the case of the mutual information. This allows to conclude that scaling 
corrections are the same for both quantities. By a standard finite size scaling analysis 
we then show that the asymptotic scaling function Rsiy) perfectly matches the CFT. 
In the Heisenberg chain both usual and unusual scaling corrections are smaller, and 
already at £ ~ 50 Monte Carlo data for Rsi^y) are in excellent agreement with the CFT 
result. 

2. Negativity and Conformal Field Theory: general results 

In the next sections we briefiy review the scaling behavior of rn{z),Rn{y) in a generic 
system described by conformal field theory (cf. [32, 33] for more details). In order to 
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make the manuscript self contained we start recalling some basic facts about Tr and 
the mutual information (which enter in the construction of in section 2.1. Then 

the behavior of r„ and Rn is discussed in section 2.2. Finally in 3 and 4 we specialize 
the result for Rn{y) to the two cases of interest: the Ising universality class and the 
free compactified boson (Luttinger liquid). The result for the Luttinger liquid has been 
derived already in [33], whereas the one for the Ising is derived here § using the results 
of Ref. [32]. 

2.1. The moments of pA (Renyi entropies) & the mutual information 

Let us consider a ID system described by a conformal field theory and take as subsystem 
A a single interval (as in Fig. 1 (a)) of length i = \vi — ui\. The asymptotic scaling 
behavior of the moments Tr of the reduced density matrix is given as 

Trpl = c„rt("-^) (11) 

with Cn a non universal constant and c the central charge of the CFT. As shown by 
Calabrese and Cardy in Ref. [3], the n-th moment of the reduced density matrix can be 
also obtained in field theory in terms of a path integral Z„ over the so-called ?7,-sheeted 
Riemann surface as 

Trp:^ = |^ (12) 

where Z is the same path integral (but on the plane) and ensures the correct 
normalization Tr p^i = 1. It is worth observing that (12) lies at the heart of all the 
algorithms for calculating Renyi entropies in both classical and quantum Monte Carlo 
simulations [15, 21, 22, 45] (as it will be better clarified in section 5). 
On the field theory side one further notices that (12) can be rewritten in terms of the 
so-called branch point twist fields T (and anti-twist T) as 

Trp:^ = (r„(ni)r.(t;i)) (13) 

The twist (and anti-twist) fields are primary fields (in the CFT language) and are inserted 
respectively at the endpoints Ui and vi of interval A (Fig. 1). Their scaling dimensions 
A„ = A„ are given as 

^n = ^{n--) (14) 



12 V 

Using (14) and basic properties of correlation functions, it is a simple exercise in CFT 
to obtain (11) from (13). 

For two disjoint intervals (see Fig. 1 (c)) the moments Trp^^^j^^ admit a similar 
representation in terms of twist fields and one now obtains the four point function 

TrPAiUAa = {%i{ui)Tn{vi)Tn{u2)Tn{v2)) (15) 

§ During the completion of this work we became aware that the same result has been derived 
by P. Calabrese et. al [44] 
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which in any CFT, using only the global conformal invariance, can be recast as 

t^^Pmua, = clliM^ - y)r^^''--^My) (16) 

Here y is the harmonic ratio (5) and J^niv) (for each n) a universal scaling function 
containing complete information about the underlying CFT. For example, the Taylor 
expansion of J-'n{y) at small y is also universal and allows to extract the full operator 
content (scaling dimensions, OPE (operator product expansion) coefficients, etc.) of the 
theory [16, 17]. From (16) the scaling behavior of the Renyi mutual information Iai-A2 
is given as 

Tr n" r n 

- log ' nt^^l = log (1 - y)-^^"-^^-^n(y) (17) 

Trp^^Trp^^ L J 

In constructing (17) the non universal factors c„ appearing in (16) cancel, implying that 
the Renyi mutual information Iai-A2 ^ universal function of solely the harmonic ratio 
y. A part from the "trivial" factor (1 — y)~'^("~i/")/^ it depends only on the function 
J^niv)- A notable consequence is that the mutual information allows to distinguish 
between CFTs with the same central charge but different operator content. The most 
prominent example is perhaps the Luttinger liquid, whose operator content changes as 
a function of the Luttinger parameter Ki, although one has c = 1 independently of Ki. 
On the other hand, one should mention that exact results for J^n{y) are known so 
far only for few CFTs, namely the free compactified boson theory [16] and the Ising 
universality class [17]. Also, even for the aforementioned models, performing the analytic 
continuation n — )■ 1 to get the von Neumann mutual information Iai:A2 represents still a 
formidable task and results are only available in some limits [16, 17]. One reason why it is 
desirable to calculate Iai:A2 is that, while IavA2 exhibit strong unusual corrections (often 
oscillating), making any attempt to extract J^niv) in microscopic models numerically 
demanding, for the von Neumann mutual information Iai:A2 scaling corrections are 
usually smaller [14, 21, 22, 23, 25, 26]. 

2.2. The moments of , logarithmic negativity, and the ratios r„ Rn 

In this section we discuss the behavior of the moments of the partially transposed 
reduced density matrix p^^^y^^ logarithmic negativity S in systems described 

by CFTs. It has been observed in Ref. [32] that Tr(p^^^y^2)"' can be written in terms of 
the same twist fields appearing in (13). In the most general case of two disjoint blocks 
Ai, A2 one has 

MptuA2r = {Tn{u{)%{v{)%MTn{v2)) (18) 

while the case of two adjacent ones can be recovered as the limit U2 — J- fi (cf. Fig 1). 
One should notice that (18) can be obtained from (15) by replacing T{u2) — > T{u2) and 
T'(f2) — > T{v2). This could be seen, somehow, as the implementation in the field theory 
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language of the partial transposition ||. From (18), using only the global conformal 
symmetry, the scaling behavior of Tr(p|J^^^J'^ can be obtained as follows [32, 33]. 

2.2.1. Adjacent intervals. We first discuss the case of two adjacent blocks, i.e. at 
distance d = (see Fig. 1 (b)). Then it can be shown that (18) reduces to a three twists 
correlation function, whose form is fully determined by the conformal symmetry. The 
result depends only on the central charge and, surprisingly, on the parity of n. Its exact 
form is given as [32, 33] 

f (44)-^^^~"^(4 + 4)-^^^+"^ neven 
MpI\uaX oc (19) 
[ (44(4 +4))-^^"-"^ nodd 

with £i(4) the size of Ai{A2) {ii given as 4 = \vi — Mil). We remark that (19) holds 
provided that 1 £±,£2 L (i.e. for two intervals embedded in an infinite chain), while 
for finite size chains one should replace, as usual in CFT, 4 L/tt sin(7r/L4) (chord 
length). The exact functional form of r„(z) now can be obtained from (19) and (6). 
Since it is quite clumsy, we do not report the explicit result, which, instead, will be 
shown numerically in Fig. 5 and Fig. 10 for respectively c = 1/2 and c = 1 (n = 3,4). 
Finally, the logarithmic negativity S for two adjacent blocks is obtained performing the 
analytic continuation n — >■ 1 of log Trdo^^)" (only using n even in (19) The result 
can be given as 

and is universal [32]. 

2.2.2. Disjoint intervals. In the (more complex) case of A being made of two disjoint 
intervals the scaling behavior of (Trp^^)" is given as 

^(pXuaJ" = c^[44(i - 2/)]-^("-^)^„(2/) (21) 

with Gn{y) a universal function of the harmonic ratio y. As for the mutual information 
(cf. previous section), the form of (21) is fixed only by the global conformal invariance 
of (18). Remarkably Gniu) can be related to the scaling function J^niy) appearing in (16) 
as [32, 33] 

^„(y) = (l-y)f("-^)j-„(y/(y-l)) (22) 

In constructing the ratio Rn{y) (cf. (7)) all the non universal factors in (21) cancel and 
one obtains a universal function of the harmonic ratio 

i?n(y) = (i-y)t("-^)^^^M^^ (23) 

II There is a subtlety here: formula (18) would correspond to Cp^C (not pj) with C the 
transformation reversing the order of rows and columns indices referring to the second interval A2. 
However, the transformation C does not affect the moments of [32, 33]. 
*\ The analytic continuation for n odd gives the normalization condition Trp^^ = 1 [32, 33]. 



Entanglement negativity and conformal field theory: a Monte Carlo study 



9 



It is interesting to investigate the asymptotic behavior of Rn{y) in the hmits y — > 0, 1. 
At y — 1 one should recover the result for two adjacent intervals: according to (5), 
in fact, y — 1 corresponds to U2 Vi. Comparing (22) with (19) one obtains that 
Tr(/?J^^^J" ~ (1 — y)'^ with 7 = c{n — l/n)/12 if n is odd and 7 = c(n/2 — 2/n)/6 
for n even. As a consequence the asymptotic behavior of Rn{y), which depends on the 
central charge and on the parity of n, is given as 

{(1 — y)T5("~n) n odd 
(24) 
(1 — y)f(t~n) n even 

On the contrary, the behavior in the limit y ^ can be obtained using the methods 
reported in Ref. [17] and is expected to depend in general on the operator content of 
the CFT. 

To conclude we mention that, since Qn{y) shows in general a non trivial dependence 
on n, it is tricky to perform the analytic continuation n — )■ 1 to obtain S. Despite 
that, according to (21) one has that, formally, £ is given as lim„_5.i log(^„(i/)) (only 
using even n, as for two adjacent blocks). This allows to conclude that the logarithmic 
negativity is a universal function of y, as already argued in Ref. [31] on the basis of 
DMRG data. Furthermore, for some CFTs exact results for £ have been worked out in 
the limit y 1 [32, 33]. A surprising result, from the CFT perspective, is that in the 
limit y ^ the logarithmic negativity apparently is vanishing faster than any power 
(i.e. in a non analytic way) [32, 33], whereas Rn{y), for any n, is analytic in y = 0. 



3. The universal ratio Rn{y) in the ID Ising universality class 



In this section we derive the universal scaling function Rn{y) in the Ising universality 
class. To this purpose we remind that for the Ising model the scaling function J-'n{y) 
appearing in (23) is given as [17] 



e(o|r 



£ 

s 



(OIF) 



where is the Riemann theta function with characteristic defined as 







£ 

L S 



fzin 



exp 



\T ■ 



(25) 



(26) 



ni[m + £j " F (m + e) + 27ri (m + e) ^ (z + 6) 

with z,£,S vectors in C"~^. Precisely, in (25) the sum is over all the possible n 
dimensional vectors £, 6 with entries 0, 1/2. The (ra — 1) x (n — 1) matrix F is defined as 



2i 



n 



k = l 



^ sin f TT- j Pk/n COS 



271— (r — s) 

n 



Here fiqix) is given by 

2Fi(g, 1 - g; 1; 1 -x) 



/3o 



2Fi(g, 1 - g; \\x) 



(27) 



(28) 
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Figure 2. Scaling function (CFT prediction) for the ID Ising universality 

class. We show Rn{y) as a function of the harmonic ratio y for n = 2, 3, 4, 5. The inset 
is to highlight the small y behavior of Rn{y). The continuous black line is ^ y^^^- 



and 2F1 is the hypergeometric function. Notice that for n = 2 (25) can be expressed in 
terms of elementary functions as [21] 

1/2 



;i + v^)(i + v^r^) 



+ xi/^ + ((l-x)x)^/^ + (l-x)i/^ 



1/2 



(29) 



The exact analytic form of Rn{y) in the Ising universality class is obtained from (25) 
and (23). This is shown numerically in Fig. 2 {Rn{y) as a function of the harmonic 
ratio y for n = 2,3,4,5). Clearly Rn{y) (for any n) is very close to one in the region 
?/ ~ and is monotonically vanishing in the limit y — t- 1. Note that R2{y) is exactly 
one (-R2 = 1 Vy). The asymptotic behavior of Rn{y) in the limit y — )■ 1, when the two 
intervals are next to each other (cf. Fig. 1), is given by (24). For instance, for n = 3,4 
one has Rsiy) ~ (1 — yY^^ and R^iy) ~ (1 — yY^^, which explains the slowly vanishing 
behavior observed in Figure 2. 

Oppositely, in the limit y — t- 0, when the two intervals are very far apart (c? ^ 1 in 
Fig. 1), the same asymptotic behavior, Rn{y) ~ 1 — for all the values of n > 2 

is observed (this is highlighted in Fig. 2: Inset). In principle the exponent 5/4 could 
be calculated analytically using the same methods employed in Ref. [17] to obtain the 
small y expansion of J^niv)- Note that one should expect a„ — in the limit n ^ 1, 
reflecting that the negativity vanishes faster than any power at ?/ — > [32, 33]. 
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4. The universal ratio _R„ in the free compactified boson theory 

In this section we re-derive (for more details cf. [33]) the asymptotic scaling function 
Rn{y) for a free compactified boson theory. This is defined by the field theory action 

S= — [ dzdzd(j)d^ (30) 



2vr . 

where are bosonic fields compactified on a circle of radius rdrde, meaning that 
= + 27ir circle- The action (30) is also the so-called Luttinger liquid, which is a 
c = 1 CFT and one of the most successful paradigms to understand the physics of 
ID systems. For instance the theory describes the critical long wavelength behavior of 
the anisotropic Heisenberg spin chain in its gapless phase, bosons with repulsive delta 
interaction, ID Hubbard model, etc.. 

For the free compactified boson theory (30) the scaling function J-'„(x) has been 
calculated in Ref. [16] and is given as 

_ e(o|ryr)e(o|r^) 
•^'^^'^^ " [0(o|r)]2 ^^^^ 

where and F are the same as for the Ising universality class (cf. previous section). 
Here r] is related to the compactification radius Vdrde as r/ = 2r^.^^;g and can be also 
given in terms of the so-called Luttinger liquid parameter Kl as 77 = 1/(2Kl). For 
n = 2 (31) is expressed in terms of the Jacobi theta functions 6^, as [14] 

where r is related to the harmonic ratio a.s y = [6'2(r)/6'3(r)]"^. Notice in both (31) (32) 
the explicit invariance under 77 — )■ l/rj. The point 77 = 4 {rdrde = V2) describes the 
critical behavior of the 2D XY model at the BKT transition (or equivalently the long 
wavelength properties of the ID Heisenberg XXZ chain at anisotropy A = —1/^2) 
(cf. [36]). 

Before proceeding one should remark that (31) is only defined for positive values of its 
argument x. Since in (23) the combination y/{y — 1) is negative for y G [0,1), one 
should consider the analytic continuation for negative argument of (31). This has been 
calculated in Ref. [33] and is given as 



n-l 



V—e(0\vG[y/{y-l)] 
J'niy) = ^ ' (33) 

n-l 

n Re(F,/„(^)F,/„(^)) 
k=i 

where Fq{x) =2 Fi{q, 1 — q,l,x) and we defined G as 
The matrix elements of G are given as 



Entanglement negativity and conformal field theory: a Monte Carlo study 12 




y Ti 

Figure 3. (a) Scaling function Rn{y) for the free compactified boson theory at 77 = 4 
{f circle = V^) ■ CFT prediction versus the harmonic ratio y for n = 2, 3, 4, 5. For each 
n the CFT resuh for the free boson in the decompactified hmit (i.e. at t] = 00) is also 
shown (rhombi). Inset: small y behavior of Rn{y), plot of 1 — Rn{y) versus y. The 
continuous black line is ^ y^^^- (b) i?„(y) (71 = 3,4,5) at fixed y = 1/2 as a function 

0f^ = 2r2^rcie- 



n— 1 |2 n—1 ^ 



k=l 



n-1 

W = — sin(7r/c/n)[sin(7r/c/n) — i cos(7r/i;/n)]i?fc/„ 

tit Nn 

where we used the definitions 

-Tfc/nia;) 

The ratio Rn{y) can now be obtained using (33) and (23). This is shown numerically in 
Fig. 3 (a) plotting Rn{y) as a function of y for n = 2, 3, 4, 5. The form of Rniy) is very- 
similar to the Ising case (Fig. 2). Precisely, Rn{y) ~ 1 in the whole region < ?/ < 1, 
while in the limit |/ — > 1 Rn{y) vanishes (faster than in the Ising case due to the larger 
value of the central charge, cf. (24)). Note, however, that Rn{y) ranges in a larger 
interval compared to the Ising case: for instance it is 0.9 < Rsiy) < 1 for < y < 3/4, 
while in the same region one has 0.98 < Rsiy) < 1 for the Ising (cf. Fig. 2). 
Remarkably, the same asymptotic behavior Rniy) ~ 1 — a'^y^^^, as in the Ising case, is 
shown in the limit ?/ — > (Fig. 3: Inset). Also it should be — )■ at n — 1 (as proven 
analytically in [33] for arbitrary values of the compactification radius). 
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Interestingly, Rniy) depends weakly on the compactification radius Tdrde, hence on rj. 
For example in the range 2 < r/ < oo, with 77 = 00 the so-called decompactified limit, 
Rn{y) does not change significantly as a function of y (in Fig. 3 the difference between 
rj = 00 and 77 = 4 is not visible at all). This is better highlighted in Fig. 3 (b) showing 
Rniy) at fixed y = 1/2 as a function of 77 = 'ir'^irde- Notice that, since Rniy) inherits the 
symmetry under 77 — )■ 1/rj from the J^niv) (cf. (31)), one has Rniv) = -Rn(l/'7) V(i/, ra). 
In Fig. 3 for each n Rniy) exhibits a minimum at 77 = 1, which corresponds to 
the antiferromagnetic Heisenberg XXX model, while in the region 77 > 1 it increases 
monotonically showing a saturating behavior for large rj. In particular already for r/ > 2 
Rniy) cannot be distinguished from its asymptotic value. The same weak dependence 
on 7] is observed at other values of y. One must mention that it is possible to calculate 
Rniy) directly in the decompactified limit (i.e. rj — )■ 00) using (33) (cf. Ref. [32] for the 
analytic result). 

5. The moments of in Monte Carlo simulations 

In this section we present a Monte Carlo scheme for calculating Tr(pJ^^^^)" using the 
replica trick and classical Monte Carlo simulations. This is based on the approach 
developed in [15, 21, 22] for Trp^ (Renyi entropies) and can be in principle applied to 
any model that can be simulated with classical Monte Carlo. The method employs the 
mapping between a quantum system in d dimensions and a classical one living in + 1. 
We should mention that an alternative numerical approach for calculating both Tr(p^^)" 
and S using TTN (tree tensor network) techniques [46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 
56, 57, 58] has become available recently [44]. 

The section is organized as follows. We first review in 5.1 the replica representation 
for the moments of both and p^^ , which lies at the heart of the method. These can 
be measured in Monte Carlo simulations using the strategy outlined in 5.2. Finally, 
in 5.3 we provide an improved scheme for models admitting a representation in terms 
of cluster variables. 

5.1. Replica trick 

The partition function Z = Tie~^^ of a d-dimensional quantum system (defined in 
terms of an Hamiltonian H) at inverse temperature /3 can be written as an Euclidean 
path integral m d+1 dimensions as 



where (^(x, r) is a field living on the hypercubic lattice {x, r} and S the Euclidean 
action. The spatial coordinates Xi are such that < Xj < Lj with i = 1,2, ... ,d and 
the imaginary time r ranges in the interval < r < = p. The fields (p are periodic 
along the imaginary time direction, i.e. 0(x, t + (3) = 0(x, r). 




(37) 
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Figure 4. Lattice version of the n— sheeted Riemann surface TZn (sec (a)) and the 
surface (see (b)) /C„ obtained sewing together n = 3 independent replicas. On each 
replica x and f stand respectively for the spatial and imaginary time directions. The 
shadow is to highlight the position of the cut, while colored links connect points on 
different replicas. We show the case of two disjoint intervals of lengths £i ~ £2 ~ 2 
at distance d = 1. In (a) (b) in each replica (plane) periodic boundary conditions on 
both sides are assumed. 



Here we consider the n— th {n G N) power of the partition function (or the rephcated 
partition function) which reads 

n " 

Z"= fU^lMe^'^'"''*'" (38) 

fc=l 

where (pk = r) is now a field hving on the k-th rephca and S((j)k) is the rephca 
Euchdean action. The actual form of the action S is not important for the following, 
but for the sake of simplicity we restrict to the case of nearest-neighbor interactions 
(which include the models treated in this work) and we consider 1 + 1 dimensions. Thus 
we assume that the action S (defined on the k—th replica) is of the form 

5(0,) = 5^F(0,(z),0,(j)) (39) 

where (ij) denotes nearest-neighbor sites and the function F models the interaction 
between the fields 0. Since we consider periodic spin chains (see Fig. 1) we assume on 
each replica periodic boundary conditions also along the spatial direction. 

5.1.1. Replica representation for the moments of Pa (Renyi entropies). We recall that 
the moments of the reduced density matrix Tr can be obtained considering the 
Euclidean partition function over the so called n-sheeted Riemann surface TZn (see [6, 4]). 
In order to define TZn we restrict for the moment to the case of A being a single interval. 
Let us start with the set of n independent replicas (sheets). Given the sheet (of area 
L X Lt) corresponding to each replica, we consider two points of its dual lattice lying 
along the spatial direction (these would correspond to the endpoints of subsystem A). 
Then we define a "cut" A as the straight line joining the two points, its length being the 
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length of subsystem A. The n— sheeted Riemann surface is defined by assuming that 
all the links starting from points on the k—th replica and intersecting the cut connect 
to points on the replica k + l(modn). 

The generalization of Tin. to multi intervals is done in the straightforward way. As an 
example in Fig. 4 (a) we show the 3— sheeted Riemann surface with L = Lr = 7 and A 
made of two disjoint intervals of equal length £ = 2 at distance d = 1. 
We now introduce the coupled action S''-"^ over the n— sheeted Riemann surface as 

n 

S^''\{M) = Y. E F{M^,Mj))+ nM^,<Pk+limodn){j)) (40) 

k=i («i>->A 
where (ij) — > A {{ij) A) denotes the links which (do not) cross the cut A. Defining as 
^n[A] = /nfc^[0fc]e"^'"'^^"^"^^ the partition function obtained from (40), finally Trp^ 
is given by [3] 

^^Pl = ^ (41) 

which is formula (12) introduced in section 2. As a final remark one should stress that, in 
order to recover Tr p\ for the one dimensional quantum system, the limit L,- = /3 — )■ oo 
has to be taken. In actual Monte Carlo simulations, however, it is sufficient to consider 
very elongated lattices with L -C 

5.1.2. Replica representation for the moments of . A representation similar to (41) 
can be obtained for the moments Tr(p^^)". To this purpose we consider a slight 
modification of action (40). Now one has A = Ai U A2, the partial transposition being 
done with respect to the degrees of freedom of subsystem A2. The cut A is made of two 
segments as A = Ai U A2, with Ai,A2 referring respectively to block Ai and A2. The 
modified action reads 



^("'^^^{0.}) = E E FiM^),MJ))+ (42) 

k=l (ij>^AiUA2 

E Fi4>k{i),(pk+limodn){j)) + E Fi^k{i),4>k-limodn)U)) 

Notice that in (42) the links crossing Ai and A2 connect fields living respectively on the 
replicas k,k + l(modrj.) and k,k — l(modn), which can be seen as the net effect of the 
partial transposition. The geometric object (that we call /C„) over which (42) is defined 
is not the ra— sheeted Riemann surface, but in general has a different topology. For the 
case n = 3 and two intervals of length £1 = £2 = 2 at distance d = 1 this is depicted in 
Fig. 4 (b). 

We checked that in our simulations for both the Ising model and XY model the choice Lr/L ^ 10 
was enough to ensure /3 = 00 within the statistical error bar (cf. also [15]) 
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After defining the partition function obtained from the action (42) as -^J^[A] one 
has [32, 33] 



Tr(p 



T2\n 



(43) 



which is the analog of (41). Both (41) and (43) can be calculated efficiently in Monte 
Carlo simulations. 



5.2. Measuring the moments of in Monte Carlo simulations 

In this section we show how to calculate in Monte Carlo simulations the ratio of partition 
functions (43) (with minor changes the same strategy applies to (41)). We start with 
defining the operator 



O = exp 



where Sx and S^^''^^'' are the "cut-linked" actions obtained by considering respectively 
in (40) (42) only the terms with (ij) crossing the cut A. Then, by definition one has 



S 



(44) 



{O) = Tr(p' 



(45) 



where (■) stands for the Monte Carlo average over the fields configurations {(f)k{x,T)}. 
These are sampled in the Monte Carlo according to the Boltzmann weights e~^'='^'''^'=^ 
constructed from the uncoupled action (39). 

At this point an important remark is in order: although the direct implementation of (44) 
in simulations is legitimate, its typical Monte Carlo history shows a huge variance, due 
to the presence of the exponential in the definition of O. This makes (44) not useful in 
practice. 



5.2.1. The increment trick. A better behaving observable, which allows to overcome 
this issue, is obtained splitting the cut A in s smaller parts A*-*^ such that A = IJ^^^^ A*^*^ 
and s is chosen arbitrarily. Defining A*-*-* = IJI=i A^^-* one has the trivial identity 



in 



11 ryT2\ 



.=0 Zi^m 

For each term in the product in (46) one can now write 

where the modified observable O is defined as 

a = exp te^^) - 5?'^^) 



(46) 



(47) 



(4^ 
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with A" = A^*"*"^^ and A' = \^^\ The Monte Carlo average in (47) is taken with respect to 
the action S^"'''^'^^ with cut A'. Now (48) receives contribution only from "fluctuations" 
living on a portion of the cut and its variance, compared to the one of (44), is strongly 
reduced. 

Further improvements are possible if the model considered (defined by the Euclidean 
action S) admits a representation in terms of clusters [59, 60] (a la Fortuin-Kasteleyn) 
and can be simulated using a Swendsen-Wang like algorithm [61]. Then it is convenient 
to express (48) in terms of cluster-related quantities. It turns out that since clusters are 
non local objects this improves dramatically the efficiency of the procedure highlighted 
so far. 

5.3. Improved Monte Carlo scheme via the "cut-linked" cluster representation 

In this section we describe the improved Monte Carlo method to simulate Tr(p^^)" (and 
Tr p^) for models that admit a representation in terms of clusters. We restrict to the 
case of 2D square lattices, although the method can be extended to models defined on 
arbitrary graphs and any dimension in a straightforward way. We first introduce the so 
called random cluster model [59]. 

Let us consider a square lattice and denote with e a generic edge (or link) connecting 
two nearest-neighbor sites x, y. Let us also define E as the set of all the links on the 
lattice and for each e ^ E consider a "link function" u{e) such that u{e) — > {0, 1}. We 
call a link e active (inactive) if u{e) = 1 (w(e) = 0), while the set of all the possible link 
configurations is Q (i.e. Q = {0, 1}^). Given an element 7 G ^2 we also consider the set 
of the active links = {e & E : uj{e) = 1}. Finally, the clusters in the configuration 7 
are defined as the connected components of C^. 

The random cluster model, which depends on the two parameters p (probability of 
activating a link) and q (cluster weight), is defined through the partition function 

Z = ^ J JJ (1 - j9)i-"(") i g^'(^) (49) 

where we denote with k{'~f) the "counting function" giving the total number of clusters 
in 7. Many models in statistical mechanics can be mapped to the random cluster 
model (Ising, Potts models, percolation models are just few well known examples). For 
instance (49) becomes the partition function of the 2D Ising model if one chooses 

q = 2 p=l-e-^ (50) 

with /5 as usual the inverse temperature. 

Clearly, the definition of the random cluster model can be extended to the surface /C„ 
(or to the n- sheeted Riemann surface TZn) in a straightforward way. To this purpose we 
first observe that the set of all the possible links configurations does not depend on the 
presence of the cut A (since the total number of links is not affected by the geometry 
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of the surface) and is given as Qjc^ = fi" = {0, l}*^^. Thus the partition function of the 
random cluster model on /C„ (i.e. Zj2[A]) reads 



g'^^^^) (51) 



where k\{'~f) is the same counting function as in (49). The subscript A is to stress that 
for a given 7 the number of clusters depends on the cut (and hence on the geometry of 
the surface), i.e. it can be kxij) 7^ kyi^) if A 7^ A'. 

In order to derive the ("improved") cluster version of (48) we observe that the partition 
function of the random cluster model on the surface /C„ with a different cut A' formally 
can be obtained from (51) (with cut A) as 



7Gf2" leG-B 



In writing (52) we used that, for each fixed 7, the term in the curly brackets does not 
depend on the cut. Dividing (52) by ^J^[A] we obtain the elegant result 

^^ = {<l-''^''h (53) 

where the average {■)\ is done (as the subscript stresses) using the action defined on the 
surface /C„ with cut A. Finally, the cluster representation of (48) is given by 

O = q''^'-^^ (54) 

with the choice A = A'-*-' and A' = A*^*"*"^^. An important property, which is useful in 
Monte Carlo simulations to speed up the evaluation of the average in (53), is that the 
difference ky — kx depends only on the clusters "intersecting" the portion of the cut 
A U A' — A n A' (i.e. only on the "cut-linked" clusters), while other contributions trivially 
cancel. 



6. Monte Carlo results: Ising universality class 

In this section we numerically investigate the properties of the scale invariant ratios 
r„ and Rn in the ID quantum Ising universality class. The Monte Carlo data that 
we present have been obtained using the results in section 5, i.e. simulating the 2D 
Ising model at the critical point /3c = l/^c = l/21og(l + -\/2). In particular we 
used the improved scheme described in 5.3 (details about the simulations can be found 
in Appendix A). 
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Figure 5. The ratios (left) and (right) versus z = t/L in the 2D classical 
critical Ising model. Monte Carlo data for several values of i = 50, 100, 150. The 
dashed dotted line is the CFT prediction (no fitting parameters). For all the data 
Monte Carlo errors are smaller than the symbols. 



6.1. Two adjacent intervals: the ratios r„ 

Let us consider two adjacent intervals Ai, A2 of equal length £ (this corresponds to 
the situation shown in Fig. 1 (b)). In Fig. 5 we show Monte Carlo data for as 
a function oi z = l/L. We considered in our simulations only L = 50,100,150 and 
1 < i < L/2. The asymptotic behavior r„(z) has been obtained analytically in CFT 
for any n in Ref. [32]. This is given in terms of (18) (cf. section 2) after replacing 
i — )■ L / 71 sin^n / Li) . The result is reported in Fig. 5 with the dashed-dotted line. 
Monte Carlo data show data collapse for both r^,r4 and for all the values of L,i 
simulated, supporting scale invariance (although this is expected only in the asymptotic 
limit L,i ^ 00). Also, the behavior of the data is perfectly reproduced by the CFT 
result. Deviations from the theory are only visible in the regions 2; ~ and z ~ 1/2. 
These are understood as finite (interval) size effects. Indeed for small z very large system 
sizes are needed to reach the asymptotic limit where CFT holds. For instance, z = 0.05 
and L = 150 (which is the largest system size we simulated) corresponds to block 
size £ = zL = 7.5, which is apparently too small and far from the scaling limit. Similar 
corrections have been observed for r^, in free bosonic systems (harmonic chain) [32, 33] 
(cf. also [44] for recent results in the Ising chain using TTN techniques). 

6.2. Two disjoint intervals: the ratio R3 

We now discuss the case of two disjoint intervals (geometry in Fig. 1 (c)) focusing on 
the behavior of the ratios Rn{y) (we restrict to n = 3). We first remind that for any n, 
in the limit L,i 00, Rn{y) are universal functions of the harmonic ratio y, and in a 
CFT are given by (22) (cf. section 3 for the analytical results in the Ising case). 
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Figure 6. Unusual scaling corrections for the ratio R-siy) in the Ising universality 
class, (top) R-siu) as a function of the harmonic ratio y (Monte Carlo data for two 
disjoint intervals of length £ — 20 and total system sizes L ~ 80, 120, 160). Different 
values of y are obtained varying the distance between the two intervals. The CFT 
prediction (dashed dotted line) is recovered in the asymptotic limit (i.e. £ — >■ oo), as 
stressed by the vertical arrow, (bottom) Plot of d3{y) = R^^^iy) — Rsiy) versus 
£-^/^. Monte Carlo data at fixed y = 0.25,0.5,0.75 and 10 < £ < 120. The dashed 
lines is the fit to da = ao£~'^^^ with og the only fitting parameter. In both plots the 
Monte Carlo statistical error bar is often smaller than the symbols. 

In Fig. 6 (top) we show Monte Carlo data for R-siy) for several values of the harmonic 
ratio < y < 1 and L = 80, 120, 160 (data at fixed i = 20). Data at different y are 
obtained varying the distance between the two intervals (according to formula (5)). 
Remarkably, all the data for different sizes L collapse on the same curve within the 
Monte Carlo statistical error, meaning that finite L scaling corrections are not visible. 
Nonetheless, at finite i = 20 Monte Carlo data do not match the theoretical (CFT) 
curve (dashed-dotted line), suggesting that finite £ corrections are present. The CFT 
result is only recovered in the limit i oo. General renormalization group arguments 
suggest the behavior 

Rn{y) = R^'^iy) + r'^"a„(z/) + • ■ • (55) 

with R^^'^{y) the asymptotic scaling function given by formula (22), while Un and 
a.„(y) are respectively the exponent and amplitude of the scaling corrections. The 
dots in (55) denote more irrelevant terms. Note in (55) the dependence of a„ on 
the harmonic ratio y. The very same behavior (upon replacing Rn{y) J^n{y) 
in (55) and fixing Un = ^/n) is shown by the scaling corrections of the mutual 
information [14, 15, 16, 18, 19, 21, 22, 24, 25, 26, 27, 17, 23]. 
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Figure 7. Amplitude 03(2/) of the unusual corrections extracted as 03(2/) = £^/^d3{y) 
with ds defined as in Fig. 6. Same Monte Carlo data as in Fig. 6 at fixed y = 
1/4,1/2,3/4. We plot 03 = i'^^^dsiy) versus ^"^Z^. The dotted lines are fits to a 
constant. Inset: fitted values for 03 (dotted lines in the main Figure) plotted versus y. 
The dashed line is a fit to Ay"^. 



More generally, entanglement related quantities are known to exhibit unusual scaling 
corrections. These are induced by the presence of the conical singularities (at the edges of 
the subsystem) needed to write the reduced density matrix (and the partially transposed 
one) in the quantum field theory language [10]. 

At difference with usual corrections, arising due to the presence of irrelevant (in the 
renormalization group sense) operators in the theory, the unusual ones are induced by a 
local insertion, near the conical singularity, of an operator that would be relevant in the 
bulk. The scaling dimension x of the operator determines the exponent of the unusual 
corrections as a;„ = 2x/n. For this reason the analysis of entanglement corrections 
provides a tool to unveil universal information about critical systems. 
On the other hand, it is a hard task in general to identify the relevant operator inducing 
the unusual corrections because any operator which does not break the symmetries of 
the model is in principle allowed [10]. For the Ising universality class it has been shown 
that this is the Majorana operator (with x = 1/2), implying Un = 1/n. It is natural 
to expect that the same one induces the scaling corrections of Rn{y), although their 
exponent could be different (i.e. Un 7^ ^/n). 

To clarify this issue in Fig. 6 (bottom) we show d^i^y) = R^^'^{y) — Rsi^y) at fixed 
y = 1/4, 1/2,3/4 and £ in the range 10 < i < 120. Clearly, Monte Carlo data support 
the behavior as We also mention that a fit to a/^'' leaving 6 as a free parameter 

gives the exponent b = 0.4(1) again consistent with b = 1/3. This allows to conclude 
that in the critical Ising spin chain the unusual corrections for Rn{y) are of the same 
form as the ones for the mutual information. 

We provide complementary information about the scaling corrections in Fig. 7, plotting 
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their amplitude a^i^y) as a function of Here a^i^y) is extracted from the Monte 

Carlo data shown in Fig. 6 as 03(2/) = £^^^d3{y). Its behavior confirms the correctness 
of the scaling corrections exponent UJ3 = 1/3. A precise estimate for as{y) is obtained 
by fitting the data with a constant (dotted lines in the Figure). The results are shown 
in the inset as a function of the harmonic ratio y and are well described by a parabolic 
behavior as ~ up to y = 3/4 (dashed line in Fig. 7: Inset). One should stress that this 
is different from what observed for the mutual information. For instance, it has been 
shown in Ref. [21] that the amplitude of the corrections to J^2{y) is very well described 
by y^^"^, pointing to a much slower decay in the limit y 0. 

7. Monte Carlo results: 2D XY model at the BKT transition 

In this section we numerically investigate the behavior of the ratios r„ and i?„ for the 
free compactified boson theory with compactification radius rdrde = (equivalently 
?7 = 4 or, in terms of the Luttinger parameter, Kl = 1/8) [36]. This describes (a part 
from logarithmic corrections) the critical properties of the 2D classical XY model at 
the BKT phase transition or, equivalently, the long wavelength behavior of the spin-^ 
quantum Heisenberg XXZ chain at anisotropy A = — 1 / -\/2 • 

The Monte Carlo data we present have been obtained using the method outlined 
in 5, simulating the 2D XY model at Pbkt (for more details about the simulations 
cf. Appendix A). Since for the XY model there is no efficient implementation of 
the Swendsen-Wang algorithm, we could not use the improved scheme provided in 
section 5.3. 

Nonetheless, one general result of this section is that the scheme described in section 5 is 
effective for models with continuous degrees of freedom. This is verified in a preliminary 
step of our analysis (cf. section 7.1) by calculating Trp^ (n = 2,3,4 with A a single 
block) and checking its scaling behavior against the well known CFT result (11). 
We then focus (section 7.2) on the ratios r3,r4 finding that, already for finite intervals, 
their behavior is well reproduced by the CFT results in Ref. [32, 33]. Surprisingly, this 
is also the case for the ratio Rsi^y) for finite (large enough) blocks size, signaling that 
corrections are smaller than in the Ising case. 

7.1. Single interval: the moments of pA CLnd the central charge 

In this section we validate the Monte Carlo procedure outlined in section 5 for models 
with continuous degrees of freedom. In particular we discuss Monte Carlo results for 
the moments Tr p\ {n = 2, 3, 4) of the reduced density matrix, demonstrating that their 
scaling behavior is fully reproduced (as expected) by the CFT result (11). Besides its 
relevance as a benchmark of the method, our analysis suggests that, even for models 
defined in terms of continuous variables, entanglement based quantities are effective 
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Figure 8. Central charge c of the Berezinskii-Kosterhtz-Thouless universahty class 
(2D XY model at the BKT point /3bkt = 1.1199(1)). We show Tip'X (Monte Carlo 
data) versus L / tt s'm^TT / Li) for n = 2,3,4, several values of L 50 < L < 200, and 
1 < £ < 75. The Monte Carlo error bar is smaller than the size of the symbols. 
Continuous lines are one parameter fits to the asymptotic CFT behavior (56) after 
fixing c = 1. Inset: Fitted values of the non universal constant c„ as a function of n 
(note the logarithmic scale on the y-axis). The dashed line is a fit to an exponential 
decay. 



tools to extract their central charge *. 

We start with reminding that the scaling behavior of Tr in the asymptotic limit is 
given in CFT as (cf. section 2) 



and c = 1 for the XY model. In Fig. 8 we show Monte Carlo data for Trp^ with 
n = 2,3,4 and 50 < L < 200, 1 < £ < 75. For each n and different L all the data 
collapse on a single curve meaning that L dependent scaling corrections are not visible 
within the Monte Carlo error bar. 

The continuous (black) lines are given by (56) where we fix the central charge c = 1, 
fitting the non universal constant c„. For n = 2 the CFT prediction reproduces the 
behavior of the MC data in the whole range 1 < £ < 75. This is not the case at n = 3, 4 
where larger i dependent scaling corrections are present and bigger systems are needed 
in order to reach the asymptotic limit. In fact, agreement with MC data is observed 

* Another available method is based on the behavior of the finite size corrections of the free energy [62] , 
which, however, is difficult to obtain accurately in Monte Carlo simulations. 
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Figure 9. The ratios rs (left) and r4 (right) versus z = i/L in the 2D XY 
model at the BKT phase transition. We show Monte Carlo data for several values of 
L = 50, 100, 150. The dashed dotted line is the CFT prediction obtained from (19) (no 
fitting parameter). For both r^^r^ the error bar is present although smaller than the 
symbols. 

only at Lc = L/7rsin(7r/L£) > 10 and > 20 for respectively n = 3,4. The same 
increasing trend in the scaling corrections of Tr ( upon increasing the Renyi index n) 
is expected in a generic system described by the Luttinger liquid [8]. 
To have an independent estimate of the central charge we performed fits leaving c as a 
free parameter in (56). For n = 2 a fit to (56) (discarding the data up to Lc < 8) gives 
C2 = 0.45(1) c = 0.99(1) with x^/DOF ^ 1.3. Similarly for n = 3, 4 one gets respectively 
C3 = 0.25(1) c = 1.00(3) and C4 = 0.13(1) c = 1.00(3). Finally, in Fig. 8 (inset) we show 
the fitted value of the non universal constant c„ as a function of n. Clearly c„ shows 
an exponential decay and is well described by the function Ae~^/'^° with A ^ 1.52 and 
no ~ 1.64. Similar exponential behavior of c„ has been already observed in the ID XXZ 
spin chain [8]. 

7.2. The ratios r„, i?„ 

We now proceed to discuss the behavior of the ratios r„. In Fig. 9 we show Monte Carlo 
data for r3,r4 versus z = IjL and system sizes L = 50, 100, 150. We considered i in 
the range 1 < £ < 75. As for the Ising universality class Monte Carlo data for different 
values of L are on the same curve (scale invariance) which is remarkably well described 
(no fitting parameters) by the CFT result (cf. (19)) (dashed-dotted line in the Figure). 
Similarly to the Ising case finite size deviations are visible in the regions z ~ 0,1/2. 
Notice also that, due to the larger value of the central charge, both r^ range in a 
larger interval, as functions of z, compared to Fig. 5. 

The behavior of the universal ratio Rsiy) is more surprising. Monte Carlo data for Rsiy) 
as a function of the harmonic ratio y are reported in Fig. 10 (data for i = 20,40,50). 
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Figure 10. Scale invariant ratio R3{y) as a function of the harmonic ratio y (Monte 
Carlo data for the 2D XY model at the BKT phase transition and i = 20, 40, 50). The 
dashed dotted line is the CFT prediction from (33). At fixed ^ = 20, L = 120 different 
values of y correspond to different distances between the two blocks. 



Focusing on the configuration with £ = 20, L = 120 Monte Carlo data are in reasonable 
agreement with the CFT curve (dashed-dotted line) for y < 1/2. Deviations from the 
theory are visible at ?/ > 1/2 and increase as a function of y, as observed in Ising model 
(cf. Fig. 6). One striking difference, however, is that apparently they decay faster 
upon increasing the block size i, which could suggest a large value of the corrections 
exponent u^. This is evident from the perfect match between theory and Monte Carlo 
for £ = 40, 50 at respectively y = 1/2 and y = 3/A. One should mention, however, that 
the precision of the data does not allow to extract a reliable estimate of the exponent 
Go's and a more careful analysis would be needed to reach a conclusion. 

8. Conclusions 

In this work we studied the moments Tr(p^^)" of the partially transposed reduced density 
matrix in critical quantum spin chains. Given a bipartition of the chains into two 
parts A and B with A being made of two equal blocks as A = Ai U A2, is obtained 
from the reduced density matrix by performing the partial transposition with respect 
to the degrees of freedom of A2. We considered both the situations with two adjacent 
and disjoint blocks, defining, respectively for the two cases, the ratios r„ and Rn{y), 
which are scale invariant at the critical point. Using a new numerical method based 
on Monte Carlo simulations we characterized their scaling behavior, in the critical Ising 
spin chain and the ID anisotropic Heisenberg XXZ model at A = — 1/a/2. 
The long wavelength properties of both models are described by well known conformal 
field theories: the first corresponds in the continuum to a free Majorana fermion (central 
charge c = 1/2), while the second is the free compactified boson at compactification 
radius rdrde = (c = !)• In li + 1 = 2 dimensions these also correspond respectively 
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to the 2D classical critical Ising model (Ising universality class) and the 2D classical 
critical XY model (Berezinskii-Kosterlitz-Thouless universality class). 
The results of our work can be summarized as follows: 

(z) Exploiting the mapping between a quantum system in d dimensions and a classical 
one in d + 1 we developed a new numerical scheme to calculate all the moments 
of the partially transposed reduced density matrix in classical Monte Carlo 
simulations. The method generalizes the one used in [15, 21, 22] for calculating the 
moments of the reduced density matrix and is effective for systems with both 
discrete (Ising model) and continuous (XY model) degrees of freedom. For models 
admitting a representation in terms of cluster variables we provided a modified 
(improved) version of the algorithm. 

{a) For two adjacent blocks Ai, A2 we studied the behavior (as a function of the length 
of the two blocks) of the ratio r„ (with n = 3, 4) in both the critical Ising quantum 
spin chain and the ID anisotropic Heisenberg XXZ mo del at A = -1/^2. In both 
cases we numerically demonstrated (for the first time in non free models fl) that 
their behavior is well described by CFT (results in Ref. [32, 33]). 

{Hi) For two disjoint intervals we studied the scaling properties of R^iy), which in the 
asymptotic limit {i — > 00) is a universal function of the harmonic ratio y. For 
the ID quantum Ising model Rsiy) exhibits strong (finite size) unusual scaling 
corrections [10]. Their exponent U3 = 1/3 is the same as that found for the 
mutual information between the two blocks [21]. After taking into account scaling 
corrections Monte Carlo data show full agreement with the asymptotic (i.e. for 
infinite blocks) CFT result. For the XXZ chain, although both usual and unusual 
corrections are expected, they appear to be much smaller and already for sizes 
£ ~ 50 we find that Rsiy) is in excellent agreement with the theory. 
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Appendix A. Monte Carlo simulations 

Appendix A. 1. The algorithm 

In this appendix we give some details about the Monte Carlo algorithms used for the 
simulations. 

jj See also [44] for a recent independent check in the Ising spin chain using TTN techniques. 
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Appendix A. 1.1. Ising model. For the simulation of the 2D Ising model we employed 
the standard implementation of the Swendsen-Wang algorithm as described in [61]. To 
generate the random numbers needed in the Monte Carlo update we used the Mersenne 
twister random number generator [63] (the same generator was used for the simulation 
of the 2D XY model). The cluster labeling step needed in the Monte Carlo update was 
performed using the standard "ants in the labyrinth" algorithm [64]. 

Appendix A. 1.2. XY model. For the 2D XY model at the BKT point our Monte Carlo 
algorithm was a local Metropolis update supplemented with some over relaxation sweeps. 
The Metropolis update that we used is different from the one usually found in the 
literature. At each Metropolis step a random vector n = (cos a, sin a) was chosen with 
a uniformly distributed in [0,27r]. Then for each spin Sx at lattice site x the update 
proposal was obtained reflecting with respect to n, i.e. 



The proposal was accepted with the standard Metropolis probability Min[e~'''^~^'\ 1] 
with E{E') the energies of the configurations before ad after the update. This procedure 
allows to avoid the (expensive) evaluation of the exponential function at each site of the 
lattice (which the is choice with the standard Metropolis implementation). 
To improve the performances of the Metropolis we added some sweeps of overrelaxation. 
The proposal for the new configuration in the overrelaxation step is obtained by 
reflecting with respect to the local field 5* = Sy (the sum is over the nearest 
neighbor sites of x): 



Since the new configuration with has by definition the same energy (as that with 
Sx) , the update step (A. 2) is accepted with probability one. However, for this reason, 
overrelaxation is not ergodic and has to be used with some other ergodic algorithm. In 
our case each Monte Carlo step was a combination of one Metropolis sweep followed by 
five overrelaxation steps. The overrelaxation sweeps have the advantage that, although 
the energy does not change, the initial and final configurations can be very different, 
reducing dramatically the critical slowing down. It can be shown indeed that the 
autocorrelation time r for the 2D XY model at the BKT point, using overrelaxation, 
can be reduced to r = O.IS.^^'^ [65] (here ^ is the correlation length) while the standard 
Metropolis gives r ~ Here it is worth stressing that any non local update (such as 
the Wolff algorithm [66]) outperforms the procedure outlined aboveff. On the other 
hand the update scheme that we used is quite effective for frustrated systems where 
there is no available cluster algorithm [67]. 

ft In particular, for the 2D XY model at the BKT point in Ref. [66], using the Wolff algorithm, no 
sign of the critical slowing down was observed. 



s'^ = 2{n-Sr,)n-s^ 



(A.l) 




2isx-S)S-s 



x 



(A.2) 
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